rm(list=ls(all=TRUE))
library(foreign)
library(reshape2)
library(plyr)
library(car)
library(ggplot2)
library(arm)
library(stargazer)
library(readxl)

options(digits=3)
as.numeric.factor <- function(x) {as.numeric(levels(x))[x]}

#get state populations from Census
census <- read.csv("state_year_population_1970_2015.csv")

#READ IN STATE ESTIMATES
sy.data <- read.csv('state_preds_means_quantiles_from_STAN.csv')
sy.data$state.op <- sy.data$mrp.trend.median
sy.data$state.op.lower <- sy.data$mrp.trend.lower
sy.data$state.op.upper <- sy.data$mrp.trend.upper
sy.data$state.year <- factor(paste(sy.data$state, sy.data$year))
with(sy.data,mean(state.op[state=="NY" & year==1993]))
n.years <-  length(unique( sy.data$year))
un.year <- unique( sy.data$year)
un.state <- unique( sy.data$state)

#now circuit data
cy.data <- read.csv('circuit_preds_means_quantiles_from_STAN.csv')
cy.data$circuit.op <- cy.data$mrp.trend.median
cy.data$circuit.op.lower <- cy.data$mrp.trend.lower
cy.data$circuit.op.upper <- cy.data$mrp.trend.upper
un.circuit <- unique(cy.data$circuit)
#cy.data.full <- read.csv('circuit_preds_all_sims_from_STAN.csv')

#now national data
nat.data <- read.csv('national_preds_means_quantiles_from_STAN.csv')
#nat.data.full <- read.csv('national_preds_all_sims_from_STAN.csv')

#policy data
policy.data <- read.csv('gay_marriage_state_bans_implementation_for_analysis.csv')

#judicial data
jud.data <- read.csv("gay_marriage_judicial_decisions.csv")
jud.data$summary <- NULL

#Which states never formally banned gay marriage?
with(policy.data, state[is.na(imp_any_ban_year)])
#CT MA NM NY RI VT
#how many states had policies set by federal courts?
with(policy.data,table(fed.court.set))

#BIG PICTURE COMBINING OPINION AND BANS/allow
states.sorted  <- arrange(subset(sy.data, year==2015), mrp.mean)$state#sort states by opinion in 2015

#add red lines for year of implementation.
min(policy.data$imp_any_ban_year, na.rm=TRUE)
table(policy.data$imp_any_ban_year)
table(policy.data$imp_any_ban_year < 1993) # 8 states
sy.data <-  join(sy.data, policy.data)#merge data
sy.data$imp_adjusted <- ifelse(sy.data$imp_any_ban_year < 1993, 1993, sy.data$imp_any_ban_year)#for graph, put pre 1993 states as of 1993

#which states never passed bans explicitly?
with(policy.data, unique(state_name[is.na(imp_stat_ban_year) & is.na(imp_con_ban_year)]))
#add Census
dim(sy.data)
sy.data <- join(sy.data, census)
dim(sy.data)
sy.data$state.pop.size.1000 <- sy.data$state.pop.size/1000
#add state-year info on whether a ban was in place
foo <- read.csv("state_year_ban_in_place.csv")
sy.data <- join(sy.data, foo)

pdf("Figure_5.pdf", height=10, width=9)

axis.text <- .9
change.text <- .8
x.adjust <- 1
y.adjust <- -.1
left.mar <- 1.5
right.mar <- .3
#par(mfrow=c(10,5),mar=c(1,2,1,1))
layout( matrix(c(1:55), nrow=11, ncol=5, byrow=TRUE) , 
	heights=c(rep(1,10),.3))
par(mar=c(.1,left.mar,.3,right.mar))
for(i in 1:50){
	ok <- sy.data$state==states.sorted[i]
	foo <- subset(sy.data,ok)
	plot(foo$year, foo$state.op, type="n", ylim=c(0,80), xlim=c(1993, 2015), axes=FALSE, xlab="", ylab="", xaxs="r", yaxs="i")
	abline(h=50, lty=2, col="black")
	polygon(c(rev(foo$year), foo$year),
			 c(rev(foo$mrp.trend.lower), foo$mrp.trend.upper),col = 'grey80', border = NA)
	points(foo$year, foo$state.op, type="l", col="black")
	move <- c("OH", "MO", "OR","NJ","MA")
 if ( ! states.sorted[i] %in% move ) mtext(states.sorted[i],3,line=-2, adj=.5)
 if (  states.sorted[i] %in% move )	 mtext(states.sorted[i],3,line=-2, adj=.2)
	#axis(1)	
	axis(2, mgp=c(2,.5,0),las=1,at=seq(0,100,25))
	box()
	#add red lines for bans
	abline(v= unique(foo$imp_adjusted), col="black", lty=2)
	#now sequentially for allow
	#legislatures or referendum
	abline(v= unique(foo$imp_allow_legref_year), col="black", lty=6)
	if(!is.na(unique(foo$imp_allow_legref_year))) text(unique(foo$imp_allow_legref_year)-x.adjust, 0, "Leg/Ref", srt=90,adj=y.adjust, cex=change.text,xpd=NA,  col="black")
	#SSC 
	abline(v= unique(foo$ssc_imp_year), col="black", lty=6)
	if(!is.na(unique(foo$ssc_imp_year))) text(unique(foo$ssc_imp_year)-x.adjust, 0, "SSC", srt=90,adj=y.adjust, cex=change.text,xpd=NA,  col="black")
	#District courts 
	abline(v= unique(foo$usdc_imp_year), col="black", lty=6)
	if(!is.na(unique(foo$usdc_imp_year))) text(unique(foo$usdc_imp_year)-x.adjust, 0, "Dist.", srt=90,adj=y.adjust, cex=change.text,xpd=NA,  col="black")
	#Courts of Appeals
	abline(v= unique(foo$uscoa_imp_year), col="black", lty=6)
	if(!is.na(unique(foo$uscoa_imp_year))) text(unique(foo$uscoa_imp_year)-x.adjust, 0, "CoA", srt=90,adj=y.adjust, cex=change.text,xpd=NA,  col="black")
	#Supreme Court
	abline(v= unique(foo$ussc_imp_year), col="black", lty=6)
	if(!is.na(unique(foo$ussc_imp_year))) text(unique(foo$ussc_imp_year)-x.adjust, 0, "USSC", srt=90,adj=y.adjust, cex=change.text,xpd=NA,  col="black")		
	}
	
par(mar=c(0,left.mar,0,right.mar))
for (i in 1:5){
	plot(foo$year, foo$state.op, type="n", xlim=c(1993, 2015), ylim=c(0,80),axes=FALSE, xlab="", xaxs="r", yaxs="i")
	axis(1, mgp=c(2,.3,0), at =seq(1993,2015,1), line=-2, labels=FALSE)
	axis(1, mgp=c(2,1,0), at =seq(1995,2015,5),las=1,line=-2, tcl=-1.2, cex.axis=axis.text)
	}
	
dev.off()

#what was opinion in MS and MA in 1993 and 2015. 
ok <- with(sy.data,state=="MS" & (year==1993 | year==2015))
with(sy.data, cbind( mrp.trend.lower[ok], state.op[ok], mrp.trend.upper[ok]))
ok <- with(sy.data,state=="MA" & (year==1993 | year==2015))
with(sy.data, cbind( mrp.trend.lower[ok], state.op[ok], mrp.trend.upper[ok]))
#VT IN 2009
ok <- with(sy.data,state=="VT" & (year==2009))
with(sy.data, cbind( mrp.trend.lower[ok], state.op[ok], mrp.trend.upper[ok]))
#hawaii in 1993
ok <- with(sy.data,state=="HI" & (year==1993))
with(sy.data, cbind( mrp.trend.lower[ok], state.op[ok], mrp.trend.upper[ok]))
#hawaii in 1993
ok <- with(sy.data,state=="HI" & (year==1993))
with(sy.data, cbind( mrp.trend.lower[ok], state.op[ok], mrp.trend.upper[ok]))

#4th circuit example
ok <- with(sy.data, circuit==4 & year==2014)
with(sy.data,data.frame(state[ok], mrp.trend.lower[ok], state.op[ok], mrp.trend.upper[ok]))

#1) create state-year data based on BAN year
policy.data$state.year <- factor(with(policy.data, paste(state, imp_any_ban_year)))
policy.data$imp_adjusted <- ifelse(policy.data$imp_any_ban_year < 1993, 1993, policy.data$imp_any_ban_year)#for graph, put pre 1993 states as of 1993
first.ban.year.data <- subset(sy.data, year==imp_adjusted)
dim(first.ban.year.data)
#timing of bans
table(first.ban.year.data$imp_any_ban_year)
table(first.ban.year.data$imp_any_ban_year<1993)
table(first.ban.year.data$imp_con_ban_year)
table(first.ban.year.data$imp_stat_ban_year)
#state op above 50 in initial implementation year?
table(first.ban.year.data$state.op < 50)
#do separately for constitutional
policy.data$state.year <- factor(with(policy.data, paste(state, imp_con_ban_year)))
con.ban.year.data <- subset(sy.data, year==imp_con_ban_year)
dim(con.ban.year.data)
#timing of bans
table(con.ban.year.data$imp_con_ban_year)
table(con.ban.year.data$state.op)

#create state-year data based on implementation year
policy.data$state.year <- factor(with(policy.data, paste(state, imp_allow_year_overall)))
imp.year.data <- subset(sy.data, year==imp_allow_year_overall)
dim(imp.year.data)

#what year did each state cross 50%?
cross.50 <- ddply(sy.data, .(state), summarise,
	year = min(year[mrp.mean>=50]))
cross.50$year <- ifelse(cross.50$year==Inf, NA, cross.50$year)
cross.50$year.impute <- ifelse(is.na(cross.50$year), 2050, cross.50$year)
median(cross.50$year.impute )
cross.50 <- arrange(cross.50, year)

#COMPARE LEGISLATIVE/INITIATIVE TO PUBLIC OPINION for legalization
	#CREATE PANEL DATA then ends until legref year
policy.panel <- subset(sy.data,select=c(state, year, state.year, state.op,state.op.lower, state.op.upper, imp_allow_legref_year, imp_allow_year_overall,state.court.set))
policy.panel <- arrange(policy.panel, state, year)
policy.panel$imp_allow_legref_year_temp <- ifelse(is.na(policy.panel$imp_allow_legref_year), 2050, policy.panel$imp_allow_legref_year)#just set to 2015 to avoid NA

policy.panel$leg.allow.gm <- ifelse(!is.na(policy.panel$imp_allow_legref_year_temp),1,0)
policy.panel$leg.allow.gm  <- ifelse( policy.panel$year < policy.panel$imp_allow_legref_year_temp , 0 , policy.panel$leg.allow.gm )
head(subset(policy.panel, state=="NM"),150)
head(subset(policy.panel, state=="WA"),150)
#now drop years after implementation
policy.panel <- subset(policy.panel, year<=imp_allow_year_overall)
head(subset(policy.panel, state=="MA"),150)
#now drop states where SSC implented
policy.panel <-  subset(policy.panel, state.court.set !=1)
#confirm it worked
table(policy.panel$state)
table(policy.panel$year)
#to drop state-years where states ultimately implemented
#ok <- with(policy.panel,!is.na(imp_allow_legref_year)& year!=imp_allow_legref_year)
#policy.panel <-  subset(policy.panel, state.court

#for non-implementing states, want policy as of year of court striking down
foo.imp.year <- subset(imp.year.data, fed.court.set==1)
foo.imp.year <- subset(foo.imp.year, select=c(state, year, state.year,state.op, state.op.lower, state.op.upper))
table(foo.imp.year$state.op>50)
#how many states have CIs that include zero
sum(foo.imp.year$state.op.lower<50 & foo.imp.year$state.op.upper>50 )
#how many above 50
sum(foo.imp.year$state.op>50)
sum(foo.imp.year$state.op.lower>50)

###########################
#MAKE FIGURE 7
###########################
pdf("Figure_7.pdf", height=5, width=8)

axis.text<-1.4
text.size <- .7
par(mfrow=c(1,1), mar=c(3,7,1,1))
sort <- -55
y.offset <- -.03
y.max <- 1
top.rug <- 1.27
plot(policy.panel$state.op,  policy.panel$leg.allow.gm,type="n",pch=19, xlim=c(30,70), ylim=c(0,y.max+.1), axes=FALSE, xlab="", ylab="", xaxs="i")
axis(1,  at=seq(10,70,10), labels=seq(10,70,10), mgp=c(2,.5,0), cex.axis=axis.text, las=1)
axis(2, mgp=c(2,.5,0), at=c(0,y.max),labels=c( "Fed. courts\nstrike ban", "Voters or\nlegislatures\nlegalize"), cex.axis=axis.text, las=1, adj=0)
mtext("State support for gay marriage", 1, line=1.5, cex=1.4)
abline(v=50, col="gray80")

# #overall logit
# model <- glm( policy.panel$leg.allow.gm~ policy.panel$state.op)
# curve (coef(model)[1] + coef(model)[2]*x, add=TRUE)

#states that legalized
	#order so can stagger labels
ok <- policy.panel$leg.allow.gm==1
xx <- subset(policy.panel,ok)

foo.vt <- subset(policy.panel, state=="VT" & year==2009) 

x <- sort(policy.panel$state.op[ok])
y <- policy.panel$leg.allow.gm[ok] [order(x)]
labels <-  policy.panel$state.year[ok]
labels <- labels[order(policy.panel$state.op[ok])]
data.frame(x,y,labels)
seq <- c(1:length(x))
ok <- seq ==1 | seq==5 | seq==8
points(x,y,type="p", pch=19)
text(x[ok],y[ok]+y.offset,labels[ok],srt=sort, adj=0, xpd=NA, cex=text.size)
ok <- seq ==2 | seq==6 | seq==9
sort<-55
text(x[ok],y[ok]-y.offset,labels[ok],srt=sort, adj=0, xpd=NA, cex=text.size)
ok <- seq ==3 
sort<- 0
text(x[ok]-.5,y[ok]-y.offset-.02,labels[ok],srt=sort, adj=1, cex=text.size)
ok <- seq==4 |  seq==10
sort<-55
text(x[ok]+.5,y[ok]-y.offset,labels[ok],srt=sort, adj=0, xpd=NA, cex=text.size)
ok <- seq==11#WA 2012
sort<-55
text(x[ok],y[ok]-y.offset-.05,labels[ok],srt=sort, adj=1, xpd=NA, cex=text.size)
ok <- seq==7 
text(x[ok],y[ok]-y.offset,labels[ok],srt=sort, adj=0, xpd=NA, cex=text.size)

text(60,.75, "Opinion in year of\nlegalization",font=3,srt=0)

#now just do opinion as of IMPLEMENTING YEAR
x <- sort(foo.imp.year$state.op)
y <- rep(.13, length(x))#set offset here
labels <- foo.imp.year$state[order(foo.imp.year$state.op)]
#points(x,y, pch=19)
#do even/odd
sort <- 0
set.seed(221)
jit.sd <- .05
ok3 <- c(1:length(x))[c(TRUE, FALSE)]
jity <- rnorm(length(x[ok3]),0,jit.sd)
text(x[ok3], y[ok3]+jity, labels[ok3], srt=sort, cex=text.size)
ok3 <- c(1:length(x))[c(FALSE, TRUE)]
jity <- rnorm(length(x[ok3]),0,jit.sd)
text(x[ok3], y[ok3]+.1+jity, labels[ok3], srt=sort, cex=text.size)
#add "legend"
text(40,.4, "Opinion in year of\nfederal court implementation",font=3,srt=10)
box()

dev.off()

#now just as 2015

#What was opinion like in:
#states where legislative passed
legref.data <- with(imp.year.data, data.frame(state=state[imp_allow_legref==1], lower.ci=mrp.trend.lower[imp_allow_legref==1], op=state.op[imp_allow_legref==1], upper.ci=mrp.trend.upper[imp_allow_legref==1]))
arrange(legref.data, op)

#states where state supreme courts struck down
ssc.data <- with(imp.year.data, data.frame(state=state[ssc_imp==1], lower.ci=mrp.trend.lower[ssc_imp==1], op=state.op[ssc_imp==1], upper.ci=mrp.trend.upper[ssc_imp==1]))
arrange(ssc.data, op)
 	#states where district courts struck down
usdc.data <- with(imp.year.data, data.frame(state=state[usdc_imp==1], lower.ci=mrp.trend.lower[usdc_imp==1], op=state.op[usdc_imp==1], upper.ci=mrp.trend.upper[usdc_imp==1]))
arrange(usdc.data, op)
#states where Courts of Appeals struck down
uscoa.data <- with(imp.year.data, data.frame(state=state[uscoa_imp==1], lower.ci=mrp.trend.lower[uscoa_imp==1], op=state.op[uscoa_imp==1], upper.ci=mrp.trend.upper[uscoa_imp==1]))
arrange(uscoa.data, op)
ussc.data <- with(imp.year.data, data.frame(state=state[ussc_imp==1], lower.ci=mrp.trend.lower[ussc_imp==1], op=state.op[ussc_imp==1], upper.ci=mrp.trend.upper[ussc_imp==1]))
arrange(ussc.data, op)

#now bring in judicial data
#note:unit of ANALYSIS IS STATE. for decision level analysis, drop duplicate citations
#for each decision, plot probability of legalization as a function of unit-specific opinion
#merge judicial data with state and circuit opinion
jud.data$year <- with(jud.data,	ifelse(cite=="191 N.W.2d 185",1993, jud.data$year))#impute 1993 opinion for Baker v. NELSOn
foo <- subset(sy.data, select=c(state, year, state.op,state.op.lower,state.op.upper))
jud.data <- join(jud.data, foo)

foo <- subset(cy.data, select=c(circuit, year, circuit.op, circuit.op.lower,circuit.op.upper))
jud.data <- join(jud.data, foo)

#create "relevant op"
jud.data$relevant.op <- with(jud.data,
	ifelse(applies == "circuit",circuit.op, ifelse(applies=="state", state.op, nat.data$mrp.trend.median[nat.data$year==2015] )) )
#now get confidence intervals
jud.data$relevant.op.lower <- with(jud.data,
	ifelse(applies == "circuit",circuit.op.lower, ifelse(applies=="state", state.op.lower, nat.data$mrp.trend.lower[nat.data$year==2015] )) )
jud.data$relevant.op.upper <- with(jud.data,
	ifelse(applies == "circuit",circuit.op.upper, ifelse(applies=="state", state.op.upper, nat.data$mrp.trend.upper[nat.data$year==2015] )) )
	
	
#impute 1993 opinion for Baker v. NELSOn
jud.data$relevant.op <- with(jud.data,	ifelse(cite=="191 N.W.2d 185", mean(sy.data$state.op[sy.data$year==1993]), jud.data$relevant.op))
jud.data$relevant.op.lower <- with(jud.data,	ifelse(cite=="191 N.W.2d 185", mean(sy.data$state.op.lower[sy.data$year==1993]), jud.data$relevant.op.lower))
jud.data$relevant.op.upper <- with(jud.data,	ifelse(cite=="191 N.W.2d 185", mean(sy.data$state.op.upper[sy.data$year==1993]), jud.data$relevant.op.upper))

#now drop duplicates 
jud.data.cut <- subset(jud.data, !duplicated(cite) ) #& (strike==fully_marriage_strike) )
jud.data.cut$state.year <- paste(jud.data.cut$state,jud.data.cut$year)
#cut.sims <- cbind(jud.data.cut,rel.op.array)

################################
#STATE COURTS--Make Figure 6
################################
foo <- subset(jud.data.cut, state_fed=="state")
#drop 1971 Minnesota case
foo <- subset(foo, cite!= "191 N.W.2d 185")
foo <- arrange(foo,relevant.op, year,state )

pdf("Figure_6.pdf", height=4, width=8)


par(mfrow=c(1,1), mar=c(3,5,.5,1))
plot(foo$relevant.op, foo$strike, type="n",pch=19, xlim=c(30,70), ylim=c(-.1,1.3), axes=FALSE, xlab="", ylab="")
axis(1,  at=seq(20,80,10), labels=seq(20,80,10), mgp=c(2,.5,0), cex.axis=axis.text, las=1)
axis(2, mgp=c(2,.5,0), at=c(0,1),labels=c( "Uphold\nBan", "Strike\nBan"), cex.axis=axis.text, las=1)
#jitter points to show confidence intervals
set.seed(2)
jitter <- rnorm(length(foo$relevant.op), 0,.1)
points(foo$relevant.op, foo$strike+jitter,pch=19)
#add 95% con. intervals
segments(foo$relevant.op.lower, foo$strike+jitter, foo$relevant.op.upper, foo$strike+jitter)
abline(v=50, col="gray80")
box()
mtext("State support for gay marriage at year of decision", 1, line=1.5, cex=1.4)

#make state symbols: by hand
#STRIKE BAN
#non-fully marriage strike
ok <- foo$state.year=="HI 1993"
text(foo$relevant.op[ok]+1, foo$strike[ok]+jitter[ok]-.1, foo$state.year[ok], srt=0, font=2)
ok <- foo$state.year=="VT 1999"
text(foo$relevant.op[ok]+1, foo$strike[ok]+jitter[ok]-.06, foo$state.year[ok], srt=0, font=2)
ok <- foo$state.year %in% c("NJ 2006")
text(foo$relevant.op[ok], foo$strike[ok]+jitter[ok]-.05, foo$state.year[ok], srt=90, font=2, adj=1)
#fully marriage
ok <- foo$state.year %in% c( "CA 2008", "CT 2008")
text(foo$relevant.op[ok], foo$strike[ok]+jitter[ok]-.05, foo$state.year[ok], srt=90, font=1, adj=1)
ok <- foo$state.year %in% c("IA 2009", "MA 2003", "NJ 2013", "NM 2013" )
text(foo$relevant.op[ok]+1, foo$strike[ok]+jitter[ok]+.08, foo$state.year[ok], srt=0, font=1)
#upholds
ok <- foo$state.year %in% c("MD 2007", "WA 2006", "NY 2006" )
text(foo$relevant.op[ok], foo$strike[ok]+jitter[ok]+.05, foo$state.year[ok], srt=-45, adj=1, font=1)

dev.off()

########
#EXTERNALITY ANALYSIS
#Use uncertainty in externality analysis. In each sim, do binomial (mean in each state, n of persons)
#EXTERNALITY ANALYSIS
#SIMULATIONS--first create "voter" database with state sizes proportional to population as of 2010
	#based on status quo before fed. courts were involved
state.data <- subset(sy.data, year==2015)
state.data <- subset(state.data, select=c(state,state.op,state.pop.size.1000))	
#state.data$op.rescaled <- (state.data$state.op-min(state.data$state.op))/ (max(state.data$state.op)-min(state.data$state.op))
state.data$phi <- state.data$state.pop.size.1000/sum(state.data$state.pop.size.1000)
state.data$floor <- 1

#need to used a different year to get states in which federal courts implemented ban
foo <- subset(sy.data, year==2005)
foo <- subset(foo, select=c(state, fed.court.set))
foo$sq.ban <- ifelse(foo$fed.court.set==1, 1,0)
foo$sq.allow <- ifelse(foo$sq.ban==0,1,0)

state.data <- join(state.data, foo)
dim(state.data)

#SIMULATIONS
#CREATE DATABASE of "voters"
set.seed(2232)#for graphing purposes
vs.dat <- state.data[rep(row.names(state.data), state.data$state.pop.size.1000/10),]#create database where "voter" is unit
dim(vs.dat)#close enough
vs.dat$state.pop.size.1000<-NULL
#set.seed(2233)
vs.dat$prefer.gm <- rbinom(nrow(vs.dat), 1, vs.dat$state.op/100)#simulate preferences

#confirm it worked
check <- ddply(vs.dat,.(state), summarise, mean.check=mean(prefer.gm))
check <- join(check, state.data)

plot(check$mean.check*100, check$state.op)
mean(vs.dat$prefer.gm)*100
mean(nat.data$mrp.trend.mean[nat.data$year==2015])
head(vs.dat,50)

###################
	#define ratio of alpha among opponents/alpha among supportes
alpha.ratio <- c(0,1)#it's a line now, so just need 2 points seq(0,1,.1)
J <- length(alpha.ratio)
total.gain.from.floor  <- rep(NA, J)

#use 1/-1 for utility

for (i in 1:J){
	#choose which type of alpha here
		vs.dat$alpha <- ifelse(	vs.dat$prefer.gm==1, alpha.ratio[i], 1-alpha.ratio[i])#as it gets higher, supporters get higher weights
		#1st do calculations without floor (See. Eq. 5 in paper)
		#in-state
			#without floor
		temp.util.in.state <- with(	vs.dat, ifelse(prefer.gm==sq.allow,1,-1))
		vs.dat$util.in.state.no.floor <- with(	vs.dat, phi*temp.util.in.state)
			#WITH FLOOR
		temp.util.in.state <- with(	vs.dat, ifelse(prefer.gm==floor,1,-1))
		vs.dat$util.in.state.with.floor <-  with(	vs.dat, phi*temp.util.in.state)
		#out-state 
			#need sum of deviations of out-of-polices from each voter's ideal point
		 #loop over states, pull out other 49 states
		 	for (s in 1:length(un.state)){
		 		ok <- vs.dat$state==un.state[s]
				state.data.temp <- subset(state.data, state!=un.state[s])
				#without floor
				state.data.temp$phi.ban <- state.data.temp$phi*state.data.temp$sq.ban
				state.data.temp$phi.allow <- state.data.temp$phi*state.data.temp$sq.allow
				vs.dat$util.out.state.no.floor[ok] <-	ifelse(vs.dat$prefer.gm[ok] ==1,
										vs.dat$alpha[ok] * sum(state.data.temp$phi.allow-state.data.temp$phi.ban), #you prefer GM 
										vs.dat$alpha[ok] * sum(state.data.temp$phi.ban-state.data.temp$phi.allow) ) #you oppose GM 
				vs.dat$total.util.no.floor[ok] <- 	vs.dat$util.in.state.no.floor[ok] +    vs.dat$util.out.state.no.floor[ok]
				#with floor
				state.data.temp$phi.floor <- state.data.temp$phi*state.data.temp$floor#redundant but clarifies
				vs.dat$util.out.state.with.floor[ok] <-	ifelse(vs.dat$prefer.gm[ok] ==1,
										vs.dat$alpha[ok] * sum(state.data.temp$phi.floor), #you prefer GM 
										vs.dat$alpha[ok] * -1*sum(state.data.temp$phi.floor)) #you oppose GM 
				vs.dat$total.util.with.floor[ok] <- 	vs.dat$util.in.state.with.floor[ok] +    vs.dat$util.out.state.with.floor[ok]
				}
			vs.dat$total.gain.from.floor <- 	vs.dat$total.util.with.floor-	vs.dat$total.util.no.floor
			total.gain.from.floor[i] <- sum(	vs.dat$total.gain.from.floor )
}

x <- alpha.ratio
plot(x,	total.gain.from.floor)


#now do over time -- loop over each year--assume 
	#a) zero alpha
	#b) constant alpha of .5 
	#c) ratio that is increasing in support over
un.years <- unique(sy.data$year)
total.gain.alpha.zero  <- rep(NA, length(un.years))
total.gain.alpha.constant  <- rep(NA, length(un.years))
total.gain.alpha.over.time  <- rep(NA, length(un.years))
sy.data.cut <- 	 subset(sy.data, select=c(state,state.op,state.pop.size.1000, fed.court.set, year, ban.in.place))	
#need to define policy in each year
mean.op <-  rep(NA, length(un.years))
mean.bans <- rep(NA, length(un.years))

foo <- arrange( subset(sy.data, year %in%c(2010,2011)), state, year)
foo <- subset(foo, select=c(state, year,state.op, ban.in.place, state.pop.size.1000))
foo$phi <- with(foo, ifelse(year==2010,  state.pop.size.1000[year==2010]/sum(state.pop.size.1000[year==2010]),
	state.pop.size.1000[year==2011]/sum(state.pop.size.1000[year==2011])))
	
for (y in 1:length(un.years)){ #1st, loop over years
	
	state.data <- subset(sy.data.cut, year==un.years[y])# & state !="CA")
	state.data$phi <- state.data$state.pop.size.1000/sum(state.data$state.pop.size.1000)
	state.data$floor <- 1
	state.data$sq.ban <- ifelse(state.data$ban.in.place==1, 1,0)
	state.data$sq.allow <- ifelse(state.data$sq.ban==0,1,0)
	mean.op[y]<- mean(state.data$state.op)
	mean.bans[y] <- mean(	state.data$sq.ban)
		
	#now create voter data
	vs.dat <- state.data[rep(row.names(state.data), state.data$state.pop.size.1000/10),]#create database where "voter" is unit
	vs.dat$state.pop.size.1000<-NULL
	vs.dat$prefer.gm <- rbinom(nrow(vs.dat), 1, vs.dat$state.op/100)#simulate preferences
	#Now alternate over alphas	
	##################################################################
		#1st: zero alpha (don't need loop)
		##################################################################
		#choose which type of alpha here
			vs.dat$alpha <- 0
			#in-state
				#without floor
			temp.util.in.state <- with(	vs.dat, ifelse(prefer.gm==sq.allow,1,-1))
			vs.dat$util.in.state.no.floor <- with(	vs.dat, phi*temp.util.in.state)
				#WITH FLOOR
			temp.util.in.state <- with(	vs.dat, ifelse(prefer.gm==floor,1,-1))
			vs.dat$util.in.state.with.floor <-  with(	vs.dat, phi*temp.util.in.state)
			#out-state 
				#need sum of deviations of out-of-polices from each voter's ideal point
			 #loop over states, pull out other 49 states
			 	for (s in 1:length(un.state)){
			 		ok <- vs.dat$state==un.state[s]
					state.data.temp <- subset(state.data, state!=un.state[s])
					#without floor
					state.data.temp$phi.ban <- state.data.temp$phi*state.data.temp$sq.ban
					state.data.temp$phi.allow <- state.data.temp$phi*state.data.temp$sq.allow
					vs.dat$util.out.state.no.floor[ok] <-	ifelse(vs.dat$prefer.gm[ok] ==1,
											vs.dat$alpha[ok] * sum(state.data.temp$phi.allow-state.data.temp$phi.ban), #you prefer GM 
											vs.dat$alpha[ok] * sum(state.data.temp$phi.ban-state.data.temp$phi.allow) ) #you oppose GM 
					vs.dat$total.util.no.floor[ok] <- 	vs.dat$util.in.state.no.floor[ok] +    vs.dat$util.out.state.no.floor[ok]
					#with floor
					state.data.temp$phi.floor <- state.data.temp$phi*state.data.temp$floor#redundant but clarifies
					vs.dat$util.out.state.with.floor[ok] <-	ifelse(vs.dat$prefer.gm[ok] ==1,
											vs.dat$alpha[ok] * sum(state.data.temp$phi.floor), #you prefer GM 
											vs.dat$alpha[ok] * -1*sum(state.data.temp$phi.floor)) #you oppose GM 
					vs.dat$total.util.with.floor[ok] <- 	vs.dat$util.in.state.with.floor[ok] +    vs.dat$util.out.state.with.floor[ok]
					vs.dat$total.gain.from.floor <- 	vs.dat$total.util.with.floor-	vs.dat$total.util.no.floor
				} #end loop over states
				total.gain.alpha.zero[y] <- sum(	vs.dat$total.gain.from.floor )
	##################################################################
		#2nd alpha=.5
		##################################################################
		#choose which type of alpha here
			vs.dat$alpha <- .5
			#in-state
				#without floor
			temp.util.in.state <- with(	vs.dat, ifelse(prefer.gm==sq.allow,1,-1))
			vs.dat$util.in.state.no.floor <- with(	vs.dat, phi*temp.util.in.state)
				#WITH FLOOR
			temp.util.in.state <- with(	vs.dat, ifelse(prefer.gm==floor,1,-1))
			vs.dat$util.in.state.with.floor <-  with(	vs.dat, phi*temp.util.in.state)
			#out-state 
				#need sum of deviations of out-of-polices from each voter's ideal point
			 #loop over states, pull out other 49 states
			 	for (s in 1:length(un.state)){
			 		ok <- vs.dat$state==un.state[s]
					state.data.temp <- subset(state.data, state!=un.state[s])
					#without floor
					state.data.temp$phi.ban <- state.data.temp$phi*state.data.temp$sq.ban
					state.data.temp$phi.allow <- state.data.temp$phi*state.data.temp$sq.allow
					vs.dat$util.out.state.no.floor[ok] <-	ifelse(vs.dat$prefer.gm[ok] ==1,
											vs.dat$alpha[ok] * sum(state.data.temp$phi.allow-state.data.temp$phi.ban), #you prefer GM 
											vs.dat$alpha[ok] * sum(state.data.temp$phi.ban-state.data.temp$phi.allow) ) #you oppose GM 
					vs.dat$total.util.no.floor[ok] <- 	vs.dat$util.in.state.no.floor[ok] +    vs.dat$util.out.state.no.floor[ok]
					#with floor
					state.data.temp$phi.floor <- state.data.temp$phi*state.data.temp$floor#redundant but clarifies
					vs.dat$util.out.state.with.floor[ok] <-	ifelse(vs.dat$prefer.gm[ok] ==1,
											vs.dat$alpha[ok] * sum(state.data.temp$phi.floor), #you prefer GM 
											vs.dat$alpha[ok] * -1*sum(state.data.temp$phi.floor)) #you oppose GM 
					vs.dat$total.util.with.floor[ok] <- 	vs.dat$util.in.state.with.floor[ok] +    vs.dat$util.out.state.with.floor[ok]
					vs.dat$total.gain.from.floor <- 	vs.dat$total.util.with.floor-	vs.dat$total.util.no.floor
				} #end loop over states
				total.gain.alpha.constant[y] <- sum(	vs.dat$total.gain.from.floor )
				##################################################################
					#3rd alpha is correlated with support over time
				##################################################################
				#use national opinion to benchmark
					nat.op <- nat.data$mrp.trend.mean[y]
					vs.dat$alpha <- ifelse(	vs.dat$prefer.gm==1, nat.op,100-nat.op)/100
			#in-state
				#without floor
			temp.util.in.state <- with(	vs.dat, ifelse(prefer.gm==sq.allow,1,-1))
			vs.dat$util.in.state.no.floor <- with(	vs.dat, phi*temp.util.in.state)
				#WITH FLOOR
			temp.util.in.state <- with(	vs.dat, ifelse(prefer.gm==floor,1,-1))
			vs.dat$util.in.state.with.floor <-  with(	vs.dat, phi*temp.util.in.state)
			#out-state 
				#need sum of deviations of out-of-polices from each voter's ideal point
			 #loop over states, pull out other 49 states
			 	for (s in 1:length(un.state)){
			 		ok <- vs.dat$state==un.state[s]
					state.data.temp <- subset(state.data, state!=un.state[s])
					#without floor
					state.data.temp$phi.ban <- state.data.temp$phi*state.data.temp$sq.ban
					state.data.temp$phi.allow <- state.data.temp$phi*state.data.temp$sq.allow
					vs.dat$util.out.state.no.floor[ok] <-	ifelse(vs.dat$prefer.gm[ok] ==1,
											vs.dat$alpha[ok] * sum(state.data.temp$phi.allow-state.data.temp$phi.ban), #you prefer GM 
											vs.dat$alpha[ok] * sum(state.data.temp$phi.ban-state.data.temp$phi.allow) ) #you oppose GM 
					vs.dat$total.util.no.floor[ok] <- 	vs.dat$util.in.state.no.floor[ok] +    vs.dat$util.out.state.no.floor[ok]
					#with floor
					state.data.temp$phi.floor <- state.data.temp$phi*state.data.temp$floor#redundant but clarifies
					vs.dat$util.out.state.with.floor[ok] <-	ifelse(vs.dat$prefer.gm[ok] ==1,
											vs.dat$alpha[ok] * sum(state.data.temp$phi.floor), #you prefer GM 
											vs.dat$alpha[ok] * -1*sum(state.data.temp$phi.floor)) #you oppose GM 
					vs.dat$total.util.with.floor[ok] <- 	vs.dat$util.in.state.with.floor[ok] +    vs.dat$util.out.state.with.floor[ok]
					vs.dat$total.gain.from.floor <- 	vs.dat$total.util.with.floor-	vs.dat$total.util.no.floor
				} #end loop over states
				total.gain.alpha.over.time[y] <- sum(	vs.dat$total.gain.from.floor )			
	}		

plot(un.years,	total.gain.alpha.zero, pch=19, type="b")
points(un.years, total.gain.alpha.constant)
axis(1, at=un.years)
abline(h=0, lty=2)

plot(un.years, total.gain.alpha.constant)
points(un.years,	total.gain.alpha.zero)			

plot(un.years, total.gain.alpha.over.time)

########################################################################
#MAKE FIGURE 8
########################################################################
#PUT ALL THREE INTO A SINGLE PLOT
pdf("Figure_8.pdf", height=13, width=7)

par(mfrow=c(3,1))

#1) POP SIZE
axis.size <- 1.4
text.size <- 1.6
title.size <- 1.2
panel.size <- 1.6

state.data <- subset(sy.data, year==2015)
state.data <- subset(state.data, select=c(state,state.op,state.pop.size.1000))	
#state.data$op.rescaled <- (state.data$state.op-min(state.data$state.op))/ (max(state.data$state.op)-min(state.data$state.op))
state.data$phi <- state.data$state.pop.size.1000/sum(state.data$state.pop.size.1000)
state.data$floor <- 1

#need to used a different year to get states in which federal courts implemented ban
foo <- subset(sy.data, year==2005)
foo <- subset(foo, select=c(state, fed.court.set))
foo$sq.ban <- ifelse(foo$fed.court.set==1, 1,0)
foo$sq.allow <- ifelse(foo$sq.ban==0,1,0)

state.data <- join(state.data, foo)
dim(state.data)

x <- state.data$state.pop.size.1000 #log(statelevel$pop.1973)

par(mar=c(3,7,2,1))
plot(x, state.data$state.op, type="n", axes=F, xlab="", ylab="", ylim=c(20,80), yaxs="i", xlim=c(0,39000))

#states that already had gay marriage
ok <- state.data$sq.allow==1 #states with bans in place before fed. courts intervened
text(x[ok], state.data$state.op[ok],state.data$state[ok], text.size, col="black", font=4)
ok <- state.data$sq.ban==1 #states with bans in place before fed. courts intervened
text(x[ok], state.data $state.op[ok],state.data $state[ok], text.size, col="black")
#lines(lowess( state.data $state.op[]~x,, f=.5), col="black", lty=2)
axis(1, mgp=c(2,.5,0), cex.axis=axis.size)
axis(2,  mgp=c(2,.5,0), cex.axis=axis.size, las=1)
mtext("State population (in thousands)",1,cex=title.size, line=2)
mtext("Estimated support for\ngay marriage (as of 2015)",2,srt=90, cex=title.size, line=2)
abline(h=50, lty=2, col="gray70")
box()
#add national opinion
abline(h=nat.data$mrp.trend.median[nat.data$year==2015],col="gray80")
text(33900, 59, "National opinion")
text(33900, 30, "A) Support versus\nstate population size", font=2, cex=panel.size)

#2 ratio
x <- alpha.ratio
y.min <- min(total.gain.from.floor)
y.max <- max(total.gain.from.floor)
plot(x,		total.gain.from.floor, ylim=c(y.min, y.max), type="l",  xaxs="i", yaxs="i", xlab="", ylab="", axes=FALSE)
axis(1, seq(0,1,.2), mgp=c(2,.5,0), labels=c(0,".2", ".4",".6", ".8",1), cex.axis=axis.size)
box()
mtext("Simulated aggregate shift in voter utility",2,srt=90, cex=title.size,line=5, font=1)
mtext(expression(paste("Ratio of sensitivity ","(",bold(alpha), ") ",  "of supporters/opponents")),1,srt=90, cex=title.size,line=2.5, font=2)
box()
abline(h=0, col="gray90",lty=1)
segments(-.1,0,0,0, col="gray90", xpd=NA, lty=2)
off <- -.02
arrows(off,0,off,y.min, xpd=NA, length=.2, lty=1, col="gray70")
mtext( "Aggregate utility declines",2, srt=90,adj=-.1, cex=.9, line=2, xpd=NA,font=3)
arrows(off,0,off,y.max, xpd=NA, length=.2, lty=1, col="gray70")
mtext( "Aggregate utility increases",2, srt=90,adj=.9, cex=.9, line=2, xpd=NA,font=3)
#label lines
#text(.56, 90839, "Supporters of gay marriage", col=support.col, srt=28)
#text(.56, 34000, "All voters", col=all.col, srt=10)
#text(.53, -44605, "Opponents of gay marriage", col=oppose.col, srt=-19)
#draw vertical
segments(.381,6.7, .381,par()$usr[3], lty=2)
text(.8, -10000, "B) Simulated utility as a\nfunction of sensitivity", font=2, cex=panel.size)
text(-.01,0,"0",xpd=NA) 

#3) over time
y.min <- min(total.gain.alpha.zero, total.gain.alpha.constant, total.gain.alpha.over.time)-100
y.max <- max(total.gain.alpha.zero, total.gain.alpha.constant, total.gain.alpha.over.time)+100
x <- un.years
text.size <-1.2
plot(x,		total.gain.alpha.zero, ylim=c(y.min, y.max), type="n",  xaxs="r", yaxs="i", xlab="", ylab="", axes=FALSE)
box()
abline(h=0, col="gray",lty=1, lwd=1)
points(x,total.gain.alpha.zero, pch=19, type="l",lty=2)
points(x,total.gain.alpha.constant, pch=19, type="l", col="black")
points(x,total.gain.alpha.over.time, pch=19, type="l", lty=6, col="black")
axis(1, at=seq(1995,2015,5), mgp=c(2,.5,0), cex.axis=axis.size)
box()
mtext("Simulated aggregate shift in voter utility",2,srt=90, cex=title.size,line=5,adj=.2, font=1)
segments(1991,0,1992,0, col="gray90", xpd=NA, lty=2)
off <- .5
arrows(par()$usr[1]-off,0,par()$usr[1]-off,y.min, xpd=NA, length=.2, lty=1, col="gray70")
mtext( "Aggregate utility declines",2, srt=90,adj=.45, cex=.9, line=1.3, xpd=NA,font=3)
arrows(par()$usr[1]-off,0,par()$usr[1]-off,y.max, xpd=NA, length=.2, lty=1, col="gray70")
mtext( "Aggregate\nutility increases",2, srt=90,adj=1.05, cex=.9, line=1.7, xpd=NA,font=3)
#label lines
text(1998, -1349, "No externalities", cex=text.size)
text(1999, -11371, "Constant externalities", cex=text.size, col="black")
text(2006, -14401, "Externalities correlated\nwith overall opinion", cex=text.size, col="black", adj=0)
text(1991.9,0,"0",xpd=NA) 

text(2012, -19000, "C) Simulated utility\nover time", font=2, cex=panel.size)

dev.off()


